#!/usr/bin/python
"""
This program makes histogram track files for Circos plots
Use MUMMER coord files as input:

Usage:
Example: dissertation_GenerateCircosTracksReadsCoverage.py ../annotation/GKIL.v6.gbf GKIL_vs_Btrimmed.1.coords GKIL_vs_Btrimmed.2.coords > ~/tools/circos-tutorials-0.56/testing/microbe.plot.illumina.txt

Note: Code is modified now to calculate both Illumina and 454 reads mummer alignments. Fix later to calculate only one if needed.

Author: Jimmy Saw
Date: 04-27-2012


"""

import sys
import re
import matplotlib.pyplot as plt
import pylab
import random
from Bio import SeqIO
from Bio.SeqUtils import GC

genome = SeqIO.read(sys.argv[1], "gb")
genome_size = len(genome.seq)

mummerfile1 = open(sys.argv[2], "rU")
m1fl = mummerfile1.readlines()

mummerfile2 = open(sys.argv[3], "rU")
m2fl = mummerfile2.readlines()

mummerfile3 = open(sys.argv[4], "rU")
m3fl = mummerfile3.readlines()

mummerfile1.close()
mummerfile2.close()

def calAvg(countlist):
    total = 0
    count = len(countlist)
    for i in countlist:
        total += i
    avg = float(total/count)
    return avg

#set initial values for genome points and adding counts
mcoords = {} #not sure why I used dictionary instead of list. maybe it's faster to set
for i in range(1, genome_size+1):
    mcoords[i] = 0

for j, line in enumerate(m1fl[4:]):
    l = line.split('\t')
    start = int(l[0])
    stop = int(l[1])
    identity = float(l[6])
    for a in range(start, stop):
        mcoords[a] = mcoords[a] + 1

for j, line in enumerate(m2fl[4:]):
    l = line.split('\t')
    start = int(l[0])
    stop = int(l[1])
    identity = float(l[6])
    for a in range(start, stop):
        mcoords[a] = mcoords[a] + 1

for j, line in enumerate(m3fl[4:]):
    l = line.split('\t')
    start = int(l[0])
    stop = int(l[1])
    identity = float(l[6])
    for a in range(start, stop):
        mcoords[a] = mcoords[a] + 1

coverage = []

for k, v in mcoords.iteritems():
    coverage.append(v)

x = 1000
depthx2 = []
depthy2 = []
while x < len(coverage):
    start = x - 1000
    stop = x
    midpoint = x + 500
    average = calAvg(coverage[start:stop])
    depthx2.append((start, stop, midpoint))
    depthy2.append(average)
    x += 1000 #increment by 1kb

cds_tracks = []
rna_tracks = []
oth_tracks = []


##For line representing the assembled genome
gx = [1,genome_size]
gy = [1,1]

xcoords = range(1, genome_size+1)

for x, y in zip(depthx2, depthy2):
    print "chr1", x[0], x[1], y

"""
##Start plotting
fig = plt.figure(1, figsize=(15,8))
ax1 = fig.add_subplot(211)

ax1.plot(gx, gy, color='#AF7817', marker='|', markersize=8.0,
    mec='black', ls='-', lw=2.0)
ax1.fill_between(depthx2, depthy2, facecolor='#33FF33', alpha=0.4)

plt.show()

"""
